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Mixed Virtual Element Methods for general second order 
elliptic problems on polygonal meshes 

L. Beirao da Veiga* * * § F. Brezzij L.D. Marini* A. Russo§ 


Abstract: In the present paper we introduce a Virtual Element Method (VEM) for the ap¬ 
proximate solution of general linear second order elliptic problems in mixed form, allowing for 
variable coefficients. We derive a theoretical convergence analysis of the method and develop a set 
of numerical tests on a benchmark problem with known solution. 

1 Introduction 

The aim of this paper is to design and analyze some aspects of the use of Virtual Element Methods 
(in short, VEM) for the approximate solution of general linear second order elliptic problems. In a 
previous paper El the same authors analyzed diffusion-convection-reaction problems with variable 
coefficients in the primal form. Here we shall deal with the mixed formulation. 

Virtual Element Methods (introduced in [TO]) belong to the family of methods that allow the 
use of general polygonal and polyhedral decompositions, that are becoming more and more popular, 
in particular in view of their use in particular problems connected to moving boundaries. Example 
of applications where polytopal meshes could have (or are already yielding) a positive impact can 
be found, for instance, in fluid-structure interaction (Ml [EE], crack propagation mmm, phase 
change SB] 05j, contact problems [22], or topology optimization 5?'i} [FT] [55], S3], but they are 
promising also in other applications, for instance in presence of coefficients that vary rapidly on 
sub-domains with complicated geometries, as when dealing with various types of inclusions (see 
e.g. [79U721 [36]), or more generally in medical applications [23 |86j [72] [77|, in image processing 
El SEES, and many others. It must be pointed out that several among these methods, in view of 
their great resistance to element distortions, come out to be handy not only for general polygonal 
elements, but also on quadrilaterals or hexahedra as well [53] . The literature on these methods has 
quite old origins (see e.g. 133), and kept slowly increasing and widening its range of applications 
ever since. See for instance si 0 ei o 0 m eh sa isa m isa eh m ma m\ ca csiisaisaiBaEii. 
In more recent times the variety of methods (already quite rich) has been growing very fast. In 
particular we have presently a flourishing group of methods, quite similar to each other, based 
(one way or another) on local polynomial reconstructions. Among others we mention Hybridizable 
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Discontinuous Galerkin methods ESI HE m mi ED, Weak Galerkin methods [Si [701 [EH ISO], the 
latest evolution of Mimetic Finite Differences mmmmm, several variants related to Finite 
Volumes and Mixed Methods [Si Hi HE1 E3 Hi Hi Hi ESj, boundary element methods m and 
various evolutions of the Virtual Element Methods themselves (mentioned below). 

The similarities and the differences among all these methods are still under investigation, as well 
as the (much more important) analysis of “which method is best suited for which class of problems”. 
We are not going to attempt to clarify these issues in the present paper, and more modestly we 
stick on Virtual Element Methods, and in particular on their use in mixed formulations. 

We recall that Mixed Virtual Element Methods for div(KV) with IK constant were introduced, 
for the two dimensional case, in [2S] as an evolution of the Mimetic Finite Differences as originally 
analyzed in [2i EB S3 Si , and then extended in various directions, see for instance 0CI3 DUE]. 
For references to several much older papers on Mimetic Finite Differences and a much more detailed 
panorama on related methods we refer to [BT] and [18] . We also point out that the first attempt to 
extend and analyze Mimetic Finite Differences to linear elliptic second order operators of the form 
div(IKV) with a variable IK was actually done earlier in [16] for the mixed formulation. 

A more recent approach to the theory of Virtual Element Methods has been introduced in [Tj, 
where the first attempt to a systematic use of the L 2 -projection operator was presented (originally 
for the so called nodal VEM). This was later refined and extended to mixed formulations in jI2] • See 
also [H] , for more details on the implementation of Virtual Elements and ns [Hi na eh EH Ezi n es 
for other interesting applications and developments. 

Here we follow this direction, and the Virtual Element Methods that we propose and analyze for 
dealing with variable coefficients are indeed based on L 2 -projection operators in a rather systematic 
way. We recall that for Virtual Element Methods the shape and trial functions are not given in 
an explicit form, but rather as solutions of PDE problems inside each element. As we do not 
want to solve these problems inside the elements (not even in an approximate way), the passage 
from constant to variable coefficients is less trivial than for other methods. In particular, simple 
minded approaches to variable coefficients can lead to a loss of optimality, especially for higher 
order methods, as it has been shown for instance in [T3] for nodal VEM. 

For the sake of simplicity we present here only the two-dimensional case, although, as pointed 
out here below in Remark |4.3] the passage from two to three dimensions, in the present case, is 
quite immediate. 

We will use the following notation. The space of polynomials of degree < k, for k nonnegative 
integer, will be denoted by Pfc, or P^O) whenever we want to stress the fact that we are working 
on a particular domain O. As common, we will use P_i = {0} as well. 

Throughout the paper, we will follow the standard notation for classical Sobolev spaces, as for 
instance in I3ZI- In particular, for a domain O in one or several dimensions, ||/||fc, p ,o (k > 0 integer 
and 1 < p < +oo) will denote the norm of the function / in the Sobolev space W k,p ((D) of functions 
that belong to L p (0) with all their derivatives up to the order k. We will also use the notation 
H k ((D) to denote W k ’ 2 ((D), and the norm of a function / in H k ((D) will be denoted by ||/||fc,e> (or 
simply ||/||fc whenever no confusion can occur). With a minor (and common) abuse of notation, 
for a vector valued function (say, f : O —> K 2 ) we will still write |jf|jfe, Pl e) to denote the norm of 
f in the Sobolev space (W k,p (0)) 2 . The scalar product in L 2 (0) or in ( L 2 (0 )) 2 will be denoted 
by (• , -)o , 0 > or simply by (• , -)o (or even (• , •)) when no confusion may arise. As usual, Hq(0) (k 
integer > 0) will denote the subset of H k {0) made of functions vanishing at the boundary dO of 
O together with all their derivatives up to the order k — 1. 

Throughout the paper, C will denote a generic constant independent of the mesh size, not 
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necessarily the same from one occurrence to the other. Sometimes, in some specific step where we 
want to stress the dependence of a constant on some variable (say, £) we will indicate it by C 
Needless to say, C £ might also assume different values from one occurrence to another. 

An outline of the paper is as follows. In Section [2] after stating the problem and its formal 
adjoint, we recall (in Subsection |2.1[ ) the mixed variational formulation. Then, in Section [3] we 
introduce the Virtual Element approximation of the mixed formulation, and derive optimal error 
estimates in Section [dj In Section [5] we derive a superconvergence result for the scalar variable, and 
finally, in Section [ 6 j we present some numerical results. 

In the bibliography we included an unusual amount of references, as it would have been appro¬ 
priate for a review paper. However we thought that a wide set of references could be convenient, 
as well, for a paper submitted for a special issue (like the present one). 


2 The problem and the adjoint problem 

Let 11 C l 2 be a bounded convex polygonal domain and let T represent the boundary of 12. We 
assume that k and 7 are smooth functions 12 —> R with /c(x) > kq > 0 for all x£ll, and that b is a 
smooth vector valued function 12 —>• M 2 . For / £ H~ 1 ( 12) (= (iLg(12))'), we consider the problem: 


J Find p £ i?d(I2) such that: 

| £p := div(—ft(x)Vp + b(x)p) + y(x)p = /(x) in 12. 

We make the following fundamental assumption, that among other things implies that problem 0 
is Well-Posed. 

Assumption WP We assume that for all source terms / £ H~ 1 ( 12) problem (|T|) has a unique 
solution p , that moreover satisfies the a-priori estimate 


lbl|i,n<C||/||-i,n, 


as well as the regularity estimate 

Iblkn < Cj|/||o,n, 

both with a constant C independent of /. 

We consider also the adjoint operator £* given by 


( 2 ) 

(3) 


£*p := div(—K(x)Vp) — b(x) • Vp + y(x)p. 


(4) 


The above assumptions on problem 0 imply, among other things, that existence and uniqueness 
hold, as well, for (92). Moreover, for every g £ L 2 ( 12) there exists a unique <p £ H 2 (12) n H^fl) 
such that £*cp = g 1 and 

IMkfi < (5) 


for a constant C* independent of g. We note that having a full diffusion tensor would not change 
the analysis in a substantial way; the choice of having a scalar diffusion coefficient ft was done just 
for simplicity. Finally, as we shall see, the 2-regularity ([ 3 ]) and ([5]) is not necessary in order to derive 
the results of the present work, and an s-regularity with s > 1 would be sufficient. Here however 
we are not interested in minimizing the regularity assumptions. 
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2.1 The mixed variational formulation 

In order to build the mixed variational formulation of problem ([Tj) , we define 

v := k~ 1 , f3 := K _1 b, 


and re-write 0 as 

u = p~ 1 (—\/p + f3p), div« + 7 p = / inti, p = 0 on F. ( 6 ) 

Introducing the spaces 

V := H( div;fl), and Q := 
the variational formulation of problem 0 is: 

{ Find (u,p) £ V x Q such that 

(vu, v) — (p, divt?) — (j3 ■ v,p) = 0 Vv £ V, (7) 

{divu,q) + ('yp,q) = (f,q) Vq £ Q. 

For the subsequent analysis it will be convenient to write ([7| also in a more compact way. For this, 
we define first 

V := V x Q, U:=(u,p), V:=(v,q), F := (0,/), 

and 

T(U, V) := (vu,v) - (p,divu) - (/3 • v,p) + (divu,?) + (7 p,q). ( 8 ) 

Problem 0 can then be equivalently written as: 

J Find Ue V such that 

|t(u,v) = (f,v) vvev. 

Remark 2.1. It is almost immediate to see that our path (from ([l])^ to can be easily reversed: 
if a pair U = (u,p) solves ([9]) then u and p satisfy ([6]) and hence p solves ([l]). In turn, this easily 
gives that the existence and uniqueness of the solution of 0 implies the existence and uniqueness 
of the solution of 0. 

3 VEM approximation 

In the present section we introduce the Virtual Element approximation of problem 0. 

3.1 The Virtual Element spaces 

Let Th be a decomposition of II into star-shaped polygons E, and let £/,. be the set of edges e of 
Th- We further assume that for every element E there exists a p E > 0 such that E is star-shaped 
with respect to every point of a disk D Pe of radius p E hE (where He is the diameter of E) and 
that the length h e of every edge e of E satisfies h e > p e He■ When considering a sequence of 
decompositions {Th}h we will obviously assume p E > po > 0 for some po independent of E and of 
the decomposition. As usual, h will denote the maximum diameter of the elements of Th- 
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For every element E we introduce: 


( 10 ) 


Gk{E) :— VPfc+i(.E), 


and 

G k (E) = the L 2 (E) orthogonal of Gk(E) in (Pfc(-E)) 2 , (11) 

so that 

(P k m 2 = Gk(E)®g£(E). (12) 

For k integer > 0 we define 

V*(E) := {v £ H( div; E) n H{ rot; E) : v ■ n\ e £ P*,(e) Ve £ dE, 

divw £ Fk(E), and rot v £ Pfc_i(.E)}. (13) 


Then we introduce the discrete spaces 

V* := {v £ R (div; fi) such that v\ E £ V/f(E) V element E in 7^}, (14) 

and 

Qh := {q £ L 2 (fl) such that: q\ E £ P k(E) V element E in %.}■ (15) 

The degrees of freedom for are obvious (one has many equivalent good choices for them), while 
the degrees of freedom for are defined by (see H21) 


f e v ■ nq k ds 

!e v ' dx 
Ie v ■aidx 


for all edge e, for all q k £ Pfc(e), 

for all element E, for all g k _i £ Gk-i(E), 

for all element E, for all g ^ £ Q^r(E), 


(16) 

(17) 

(18) 


where the notation (pf0[)-( 111 was used for G k (E) and G^(E), respectively. 


Remark 3.1. We point out that conditions (161 could be replaced by the values of v ■ n at suitable 
points on each edge. Similarly, in (18) G k (E) could be replaced by any subspace of (P k (E)) 2 
satisfying (12). 


Remark 3.2. It is not difficult to check that the present choice of elements mimics, in some sense, 
the Raviart- Thomas elements, although, even on triangles, they coincide with the RT elements only 
for k = 0. j4s pointed out in lTf9J and in \12j there are many other choices that could be made. 


Remark 3.3. Regarding the mesh assumptions at the beginning of this section, we note that it 
wouldn’t be a problem to generalize the shape regularity condition by allowing suitable unions of 
star-shaped elements. Analogously, also the minimal edge length condition could be probably avoided 
with some additional technical work in the interpolation estimates. 


3.2 Interpolants, projections and approximation errors 

From now on, we shall denote by 11^ : Q —»• Q'f and by 11° : V —> Vjf the L 2 — projection operators, 
defined locally by 


/ (q ~ n ° k q)p k dx = 0 \/p k £ P k (E), ME £ %, 

J E 

[ (v - n ° k v)q k dx = 0 Mq k £ (P fc (.E)) 2 , ME £ T h . 
J E 


(19) 
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In P2] it was shown that the degrees of freedom (16)-(18) allow the explicit computation of the 
projection 11 °v from the knowledge of the degrees of freedom (16)-(18) of v. For the convenience 
of the reader we briefly recall the construction. We first observe that using the degrees of freedom 
(17) we can easily compute the value of v ■ n on dE. From this and (17) one can compute the value 
of div v £ Pfc, using 


div v q k dx = — / v ■ 'S/qkdx + / v-nq k ds \/q k £ 


( 20 ) 


I dE 


(remember that Vgt € Gk- i). Once you know explicitly v ■ n on dE and div v inside E, then you 
can easily compute the integral 


v-\7q k+1 dx=- / div v q k+ \dx + / v-nq k+1 ds, 


IE 


( 21 ) 


meaning that you can compute f F v ■ g k dx for every g k £ Gk- This and the degrees of freedom (18) 
allow you to compute f E v ■ q k dx for every (vector valued) polynomial q k of degree < k. 

On the other hand, in every element E, the computation of the L 2 (i?)-projection of an element 
q £ Q k is trivial (and coincides with its restriction to the element E). 

With classical arguments one can easily show that 


\q-n° k q\\o<Ch s \q\ s , || V - H°v\\ 


< Ch s \v\ s , 0 < s < k + 1, 


( 22 ) 


for every q and v , respectively, that make the norms in the right-hand sides finite. 

We point out that a linear “Fortin” operator 11^ from W := (H ^H)) 2 
through the degrees of freedom (16)-(18), by setting, brutally 


V k can be defined 


I e ( v ~ u h v ) -nq k ds = 0 
J E ( v - U h v )-9k-idx = Q 
J E ( v - U h v )-9i dx = 0 


for all edge e, for all q k £ Pfc(e), 

for all element E, for all g k _ 1 £ G k -i(E), 

for all element E, for all g k £ G k (E), 


(23) 

(24) 

(25) 


and (using, essentially, (20)) it is easy to verify that the commuting diagram property holds: 

div 


n 


w 

f) 

vt 


-* Q 


n° 


div 


-* Q 


0 


0 


so that 


div Tl k v = 11° div v. 


(26) 


(27) 


Moreover, the following estimates hold, provided u has enough regularity: 

||u - n£u||o < C'/i fc+1 ||M|| fe+ i, || div(w - n^tt )|| 0 < Ch k+1 \\ divtt|| fc+ i. (28) 

With a minor abuse of notation, for an element W = (w,r) with w £ (Ffg(f2)) 2 and r scalar or 
vector function in L 2 (fl), we will also denote 


w := n ° k w, 
Wi := Tl^w, 


r := II°r, and W := (w,r), 
ri := n°r, and Wj := (tt>/,r/). 


(29) 
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We remind that, obviously, 


w\\o < Kilo, 


nlo < IMIo, 


Vi llo < 


(30) 


while 


K/llo < Iloilo + IK - ™/||o < C(|K||o + /iHi). 


( 31 ) 


For locally smooth w, as we can see from (22) and (28), the two errors \\w — tu/||o,B and 
|K — m||o,_E will behave in the same way (in terms of powers of h and required regularity). Hence 
it makes sense to introduce a sort of common value that bounds both of them. We define 


£ k {w) := K - ™/||o + IK - Hlo, 


(32) 


and we put it on charge to measure the approximation error for w when using Virtual Element 
spaces of degree k. Needless to say, the same holds for a scalar function r approximated in Q k (or, 
when necessary, for a vector valued function r approximated in ( QV) 2 ) since, in these cases the two 
approximations r and ry coincide (as one can see in (29)). In order to use the same notation all 


over, however, we follow (32), and set 


\r) := llr-r/Ho + llr-rllo, and K(W) := ||W - W 7 || 0 + ||W - W| 


(33) 


We also point out that, by the properties of the projection, we immediately have 

||W/ - Wj||o < ||W/ - W||o < ||W 7 - W|| 0 + ||W - W|| 0 = £ k (W), (34) 


implying also 

£ k (Wi) < K(W). (35) 

Along the same lines, it is intuitively obvious (and it can be easily proved) that if you have a 
certain estimate (in terms of powers of h and required regularity) for w (or for r) you will have 
quite similar estimates for, say, ipw whenever ip is a given smooth function. The constant in front of 
the estimate will depend on ip, but the power of h and the regularity required to w will be exactly 
the same. For instance it is immediate to check (just expanding the derivatives of the products, 
and using Cauchy-Scliwarz) that one has 


\\pw - vKlo < Ch k+1 \ipw\k+i < Ch k+1 |MU+i,oo|Kllfc+i = IKIU+ 1 - ( 36 ) 


The same occurs for a pair W = ( w,r ) when one of the two entries (or both) are multipled by 
a smooth function <p or a smooth vector valued function ip, as in 

£ k (wtp) = \\wp - wtp\\ Q + \\wp - KvO/Ho and £k ( ri P) = Wv-rvh + Wv- (^)/||o (37) 


for a smooth function ip, as well as in 

£ k (rp) = ||rp - np\\ 0 + ||rip - (rv?)j||o (38) 

for a smooth vector valued function ip. 

All this suggests a further “abuse of notation”: for W = (w, r) we will use the notation f k («W) 
(either for n scalar or « vector) whenever one of the two (w and r), or both, are multiplied by n. 
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It could be worth pointing out a few particular cases: no matter whether « is a scalar or a 
vector, we have 

(39) 


^("W) < C^h ^||r|| ! + H^lli), 


as well as 


£ fe (xW/) < £ k («(W! - W)) + £' fc (xW) < |M|oo£ fc (W) + £ k {«W). 


Needless to say, the obvious analog of the bounds 
£ k (w), £ k [r ) and so on. Finally, we observe that estimates (22) and (28) imply 


£ fc (U) < Ch k+1 (\\u\\ k+1 + ||p|| fc+1 ), £ k («V) < C«h k+1 ( |M| fc+1 + |b|| fe+1 ), 


(40) 

-(40) apply also to the separate terms 

(41) 


where C N is a constant depending on n and its derivatives up to the order k + 1. As a final remark 
we note that, whenever convenient , we can easily bound £ k (£W) by 

£ k (nW) < ||W|| 0 - (42) 


3.3 The discrete bilinear forms 


As is well known from the theory of mixed formulations, the two main ingredients to be used to 
prove stability and error estimates are the ellipticity of the leading diagonal term (here, (vu,v)), 
and the inf-sup condition. Here the inf-sup condition will be easily provided by the commuting 


diagram (26). Hence, our main worry will be the treatment of the term 

a{u , v) := {vu 1 v). 


(43) 


On each element E £ Th we define: 

a E (v, w) := {vv, w)o,e + S E (v - v, w - w), (44) 

where S E (v,w ) is any symmetric and positive definite bilinear form that scales like a E {v 1 w) (see 
m). More precisely, our assumption on S will be: There exist two positive constants a* and a* 
(depending on v but independent of h ) such that 

a*a E (v 1 v) < S E (v, v) < a*a E (v 1 v) \/v&V k . (45) 


For practical purposes it will be convenient to choose the Euclidean scalar product associated to 
the degrees of freedom in V k multiplied, for instance, by |A|^(xs), where xb = ( Xb,Vb ) = is the 
barycenter of E. We notice that, obviously, p fc = p& for all Pfc € Pfe. Therefore 

a E (p kl w) = I j/p fc w, dx Vw £V k , Vp k £P k - (46) 

J E 

We can now define 

a h {v,w) := J2 a h(v,w). (47) 

E 

Lemma 3.4. The bilinear form ah{-,-) is continuous and elliptic in (L 2 (H)) 2 , that is: 


3M > 0 such that \ah(v, m)| < M||i>|| 0 |M|o Vv,weVfi, 
3a > 0 such that ah(v,v) > ck||||§ € V k , 


(48) 


with M and a depending on v but independent of h. 








(49) 


Proof. The symmetry of S E and ( |45| imply easily the continuity of S E : 

S E {v, w) < (S E (v, v)) 1/2 (S E (w,w)) 1/2 < C i/ |M| 0 ,.e|M|o,.e, 
with C v = a*v max . In particular, 

S E (v -v,w-w)< C v \\v - v|| 0 ,£)||to - «J||o,b < C v £ k (v)£ k (w). 


(50) 


Then, the continuity of dh(-, •) is an obvious consequence of the continuity of a(-, •) and of the L 2 — 
projection properties: 


\ah(v,w )I < twIMIolMlo + C v \\v - v\\ 0 \\w - tu|| 0 < M ||u|| 0 ||uj|| 0 . 


Similarly, 


ah(v,v) > t'min(||v || 2 + a*\\v - u||o) > a(||u 


= r " 2 + ||«-tJ||^) = allvlln- 


The discrete problem is now: 


{ Find (■ Uh,Ph ) £ V k x Q k such that 
a h {u h , v h ) - ( p hl divv h ) - (/3 • Vh,p h ) = 0 Mv h e V k 
(div u h , q h ) + ( 7 p h , q h ) = (/, Qh) Vqh, £ Qh- 


□ 


(51) 


Like we did for the continuous formulation, in order to write (|51j) in a more compact form, we set 
V h :=V k xQ k , U h ~(u h ,p h ), V/j := (y, q h ), F„ := (0 ,/), 

and 

A(U h ,V h ) := a h (u h ,v h ) - (p h , div v h ) - (/3 ■ v h ,p h ) + (di vu h ,q h ) + (7 Ph,Qh)- (52) 

Then problem (|51|) can be written as 


Find Vh £ Vh such that 

A h (V h ,V h ) = (F hl V h ) VV fc e V h . 


(53) 


4 Error Estimates 


Our final target is to prove the following theorem. 

Theorem 4.1. Under the above assumptions and with the above notation, for h sufficiently small 
problem (51) has a unique solution ( Uh,Ph ) £ V k x Q k , and the following error estimates hold: 

\\p-p h \\o<Ch k+ 1 (\\u\\ k +1 + IblU+i), 

\\u-u h \\o<Ch k+1 (\\u\\ fc +1 + Ibllfe+i), (54) 

|| div(w - u h )\\ 0 < Ch k+1 [\f\ fc+l + INI fc+i) > 

with C a constant depending on v,f3, and 7 but independent of h. 

Before proving the theorem, we will introduce some useful lemmata, that deal with properties 
of the bilinear forms A and Ah- 
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4.1 Preliminary estimates 

A typical source of difficulties, when proving optimal error estimates, is the fact that the bilinear 
form _4(U,V) cannot be bounded in terms of the L 2 norms of U and V, due to the presence 
of the two terms (di vu,q) and (p, divu) involving the divergence. We will therefore spend some 
additional time in order to point out some particular cases in which these terms could be avoided. 
In particular, we note that for v £ H{div;E) and q £ L 2 {E) we will have 


div v qdx = 0 


JE 

whenever 

• q £ Pfc, and divv is orthogonal to P&, 

• div « £ Pt, and q is orthogonal to P^. 


Hence, in particular, using (13), (15), and (27) we have: 


/ div(w - Yl^w)q h dx = 0 \/q h £ Q k h , V w £ {H 1 {E )) 2 , 
J E 


and 


f div-^ (r - n° k r)dx = 0 \/v h £ V£(E), Vr € L 2 {E), 

JE 

so that for every W £ V and for every V/, £ V;, we have 


(55) 


(56) 

(57) 


|^(V h ,W-W / )| + |^(W-W / ,V h )|< a i /9 l7 ||V h ||o||W-W / || 0 . (58) 


4.2 The consistency error 


Further attention should also be given to the difference {Ah — -4)(W, V). We will perform the 
analysis on a single element, without indicating every time that the norms are considered in L 2 {E). 
Using (|8]) and (52) we have easily 


{Ah - A)(W,V) = (vw, v) - {ew,v) (=: Ti(W, V)) 

+ S{w-w,v- v) {=: T 2 (W, V)) (59) 

+ {v-v,Pr) (=: T 3 (W, V)), 


where as before V = {v, q) and W = (w,r) are in V h - We point out that all the terms T 1; T 2 and 
T 3 do not involve derivatives, so that we will not have continuity problems. For the term Ti, using 
repeatedly the properties of the L 2 — projection we have: 

Ti(W, V) = ( vw , v ) — {ew, v) = {ew, v — v) — {w — w, ev) 

= (EW — EW 1 V — v) — {W — W, EV — EV ) 

= {eW — EW , V — v) — {w — W, EV — EV + EV — Ev) 

= {ew — ew 1 v — v) — {w — w, ev — ev) — {w — w , e{v — u)) (60) 

< (c„\\w - uJ||o + II EW - FuT||o) ||u||o 

< (a£ fc (W)+f fc (^W))||V|| 0 . 
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Needless to say, in view of the symmetry of the term, we also have 


T 


(W,V) =T 1 (V,W) < (c v £ k (\) + £ k (v\))\ 


W|| 


(61) 


The terms T 2 and T 3 in (59) are easily bounded. Directly from (|50j) we have 

T 2 (W,V) < C v £ k (W)£ k (V), 

and for T 3 

T 3 (W, V) = (v- V,(3r-pr) < £ k (V) £ k (/3W). 


(62) 


(63) 


The above analysis can now be summarized in the following two estimates, that will both be used 
in our final proof. 


Using (61), (63) with (42), and (62) we have 


(A h - -4)(W, V) < C v , p ((£ fc (V) + £ fc (zzV)) ||W|| 0 . 


(64) 


Using instead (60), (63) and (62) we have 


(A h -A)( W,V) < C v (£ fc (W) + £ k (vW) + £ k (/3W)^j ||V|| 0 . 


(65) 


4.3 The dual problem 

Our proof will use a duality argument. Therefore we spend some time analyzing the dual problem. 
Lemma 4.2. Let £ £ L 2 (£l), g £ H( div; £1), and set G := ( g , £). Let Z := (£, z) £ V be the solution 


of 

A( W,Z) = (G,W) V W = (w, r) £ V. (66) 

Then Z is the solution of 

C = K(Vz + g) and — div (, — (3 ■ £ + 7 z = £ in £ 1 , z = 0 onT (67) 

that is, (see 

2,*z = £ + b • g + div ( «g), ( 68 ) 

so that, in particular 

INI 2 + HCII 1 <C*(||4|o + ||Kfl||ff(div)). (69) 

Proof. Recalling (| 8 ]), and substituting W for U and Z for V we get 

*4(W,Z) = (vw,Q - (r,divC) - (/3 • C,r) + (div w,z) + (7 r,z). (70) 


Separating the equations in w and in r in ( |66| ) it is not difficult to see that (£,z) solves 

f (vw, Q + (div w,z) = (g, w) \/w £ H (div, £1) 

[ - (r, div 0 - (/3 ■ C> r) + ( 7 r, z) = (£, r) Vr £ L 2 (£l) 
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giving, respectively, 


C = kVj + Kg plus z £ Hq(SY), 


and 

— div C — /3 • C + l z = t- 

Putting them together we have 


— div(/cV:r) — div(Kg) — b • V z — b ■ g + 'yz = £, 


and ( 68 ) follows. 


□ 


4.4 Proof of Theorem 14.11 

We are now ready for the proof of Theorem |4.1 


Proof. To prove Theorem 4.1 we shall follow the arguments of Douglas-Roberts m ■ We first 
assume that (51) has a solution, at least for h sufficiently small. That it does, it will be clear from 
the convergence analysis. Let therefore XJh = (uh,Ph) be a solution of (51). Let us form the error 
equation: 

A(U,V h )-A h (U h ,V h ) = 0 VV* = K ,?n)6Vi. 


We use duality arguments. Let M/ = (x, if) be the solution of the adjoint problem 
A(V,*)= (i/(U/-U h ),v) = -u h ), Pl -p h ),v) We V. 


According to Lemma 4.2 if e D iJ 2 (fl) is the solution of the adjoint problem 

£*'!/’ = div(-WV’) - b- V'ip + r yip=pi-ph + /3 - (uj - u h ) + div(it/ - u h ), 
and by the elliptic regularity (69) with G = (g,H) := (v(uj — Uh),pi — Ph) we get 
HV’lb + ||x||l < C* (|| P! - Ph||o + || Ul ~ U h |\ H (div))- 


(72) 

(73) 

(74) 

(75) 


Our first step will then be the estimate of || div(tij — Mh)||o- Looking at the discrete and continuous 
equations we have 


div u h = U° k (f - 7 p h ) and divu = / - yp, 


and from (27) div it/ = 11° divit = 11° (/ — yp). Hence, 


so that, clearly, 


div (it/ - u h ) = n°(y(p/j - p)), 

|| div(U/ -Mh )|| 0 < C'-y||p — Pft, ||o- 


Therefore, (75) reduces to 

Hh + llxlli < C(||P/ — Philo + \\ui - Uh||o + ||p — Pz||o) 

< c(l|U/ — Uh||o + £ fe (U)), 


(76) 

(77) 

(78) 

(79) 
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and using (for instance) (39) with n = 1 the estimate (79) implies that 

£*(¥)< Cfc(||V>||i + ||x||i) <Ch (||U/-U h ||o + £ fe (U)), 


as well as 


\*i\\o < 11* - */||o + Iloilo < c(||U 7 - V h \\o + £ k (U))- 


Moreover, taking V = U 7 — in (73), it is immediate to see that 


Hence, 


|U 7 -U h || 2 < f (v\uj - u h \ 2 + \ PI - Ph \ 2 )dx = A(U T -U h) ¥). 
■hi 

^min||U/ - U J 2 < M(U 7 - U„, *) (±* 7 ) 

= M(U 7 - U h , * - * 7 ) + A(Ui - U h) * 7 ) (±U) 

= I + _4(U 7 — U, \I> 7 ) + A (U — U/j,\& 7 ) (just linearity) 

= I + II + M(U,* 7 )-M(U ft ,* 7 ) (use@) 

= I + II + {Ah — ^4)(U/j, ’3'/). 


( 80 ) 

(81) 

(82) 


(83) 


The first two terms are easily bounded using (58), (80)-(811, and (41): 

I = A( U 7 - V h , * - * 7 ) < C ||U 7 - U fc || 0 /i(||U 7 - U fc ||o + £‘ fe (U) 

< c(h\\Ui - U^ll 2 + h k+2 \\XJi - Ufcllo), 

II = M(U 7 - U, * 7 ) < C£ k ( U) (||U 7 - Uhllo + £ f '’(U)) 

<c(||U 7 -U,J| 0 /i fe+1 + /i 2fe+2 ), 

and we are left with the third term. For it, we are going to use the arguments of Subsection 4.2 
We start by observing that 


(84) 


(85) 


(A h - A)(\J h , * 7 ) = (A h - A)(U h - U 7 , * 7 ) + (A h - -4)(U 7 , * 7 ). 


( 86 ) 


The first term in (|86j) can be easily bounded, using ([64]), (|35j), (|40j), (|80j), and ( |41[ ): 

(A h - A)(U h - U 7 , ® 7 ) < C Vll3 (£ k (^i) + £> V*/)) ||Uh - U 7 1|o 

<Ch (||U^ - U 7 || 0 + £ fc (U)) ||U h -U 7 || c 
<c(h ||U fc - U 7 ||q + h k+2 \\\J h - U 7 1| 0 ), 


(87) 


while, using ( |65| ), (351, (40), and (81), the second term in ( |86[ ) can be bounded by 
(A h - A)(Ui, * 7 ) < C v (f fc (U 7 ) + £ k (u U 7 ) + £ fc (/3U/)) ||* 7 || 0 

C (>(U) + £ k (uU)) + £ fe (/3U)) (||U h - U 7 1|o + £ k (U) 
c(/i fe+1 ||U fe -U 7 || 0 + /i 2fe+2 ). 


< 

< 


( 88 ) 
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Inserting (84), (85), and |87|)-(88) into (83) we have then 

Z'minllU h - U/ll 2 < C (/l||Ufc - U 7 f + ||Ufc - U/ll h k+1 + h 2k+2 ) . 

For h small enough (say: Ch < (1/2 )v m in in ©) we can hide the first term in the r.h.s. of 
in the left-hand side, and have 


( 89 ) 


IU^-UjII 2 ^ c(ft* +1 ||U fc -U/|| 0 +/i 2fe+2 ), 


(90) 


and the first two estimates in (54) follow completing the square. The estimate on the divergence 


follows directly from (76), and standard error estimates. 


Finally, since (51) is finite dimensional, in order to prove the existence of the solution we only 
have to prove uniqueness, that is, we have to prove that for f = 0 problem (51) has only the solution 
Ph = 0, Uh = 0. Since we assumed that the continuous problem (|T]) has a unique solution, it follows 
that for / = 0 we have p = 0,u = 0. The above analysis showed that, for h small enough, any 


solution (v,h,Ph) °f (51) must satisfy (54) which, in our case, imply Uh = 0,ph = 0, and the proof 
is concluded. 

□ 


Remark 4.3. Looking at the construction of the method, and to the analysis of its convergence 
properties, it is not difficult to see that the passage from the two-dimensional case to the three- 
dimensional one can be done, using m without any difficulty. However, the notation for dealing 
with both cases at the same time would be more cumbersome, and a presentation with two separate 
treatments would be very boring and essentially useless. 


5 Superconvergence results 


Theorem 5.1. Let ph be the solution of (51), and let pj G Q k be the interpolant of p. Then, for 
h sufficiently small, 


\\Pi~Ph\\o < Ch k+2 (\\u\\ k +1 + HpIU+i + l/U+i 


(91) 


where C is a constant depending on v , /3, and 7 but independent of h. 

Proof. We proceed again via duality argument. Let ip € D H 2 (Ll) be the solution of the 

adjoint problem 


div(— k(x)V^)— b(x) • V^+7(x) ip = p 7 - p h , % = nVip, 
whose mixed formulation is: Find (x,V0 i n H{ div, fl) x L 2 (fl) such that 

f v) + (ip, div v) = 0 \/v € H (div, f l) 

I - (divx,g)-(/3 ■ X,9)+(7V’,9) = (pi ~Ph,q) Vg G L 2 (Li). 


The error equations 


using (29) and (27), become 


a(u,v h ) - a h (u h ,v h ) - (p 7 - p h ,divv h )-((3 ■ v h ,p) + ((3 ■ v h ,p h ) = 0 Mv h G V k , 
(div(w - u h ),q h )+fy(p-p h ),q h ) = 0 Vq h G Q h . 


(92) 


(93) 


(94) 
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Taking now q = pi — Ph in (93) gives 


\\Pi~PhWo = -(divx,_p/ - p h )-((3 ■ X,Pi ~ Ph)+hi>,Pi ~ Ph)- 


(95) 


For the first term, using again (29) and (|27j) , we have 


(div XiPi ~Ph) = (divx/,p/ - Ph) (use(|94]) with v h = Xi) 

= a ( u i Xi) - a h (u h ,Xi)-{P ■ Xi,P) + {P ' XiiPh) (±u h ) 

= a(u - u h , Xi) + a(u h , Xi) - a h (u h , Xi)~(P • Xi,p) + (P ' Xj,Ph) 
= a(u - u h , x) + a{u - u h ,Xi - x) + a(u hl Xi) - a h {u h , Xi) 

~(P ■ Xi,p) + (P -xi,Ph)- 


(±X) (96) 


In turn, the first term in (96) becomes 

a(u - u h ,x) = (u- u h ,Vip) = — (div(it - u h ), 'ip) (±ipi) 

= —(div(it - u h ), ip-ipi)- (div(it - u h ), ipi) (use ( f9lj )) 
= ~(div(u-u h ),ip -ip^+iyip^p-ph)- 


(97) 


Replacing (971 in (96), and using the result for the first term of (95), we have then 

\\pi - Philo = - [ - (div (it - u h ), ip - ipi)+{'yipi,p - p h ) + a{u - u h , Xi ~ X) 
+ a(u h , Xi) - a h (u h , Xi)~{P ■ Xi,p) + {P ' Xi,Ph)\ 

~(P ■ X,Pi -Ph)+(r>p,pi -Ph)- 

The first two terms are easily bounded: 

I a(u - u h ,xi - x)l < Ch\\u - ithllolb/ -Philo, 

|(div(u — Uh),ip — ipi)\ < Ch 2 \\ div(u — ith)||o||p/ — Ph||o, 


(98) 


while using (611 and ( |62[ ) we get 

I a(u h ,Xi) - a h (u h ,xi )I < C v h k+1 \\u\\ k+i,nh ||pi — Ph\\o- 


(99) 


( 100 ) 


For the terms involving reaction, adding and subtracting (yipiiPi — Ph) and using the properties of 
the projection we obtain 


(')ip,pi-ph)-(jipi,p-p h ) = (l(ip - ipi),Pi ~Ph) + (7 i>i,Pi~Ph) - (t>Pi ,P~Ph) 

= (t(V ) - i>i),Pi ~Ph) + (7 i>i,Pi~P) 

= ( 7 (ip - 4>i),pi - Ph) + (-yip 1 - yip 1 , pi -p) 

< C^h 2 (\\p! -Ph||o + ||p - P/||||P/ -Ph||o)- 


(101) 


For h small enough the first term in the right-hand side of (101) can be hidden in the left-hand 
side of (98) and the other one is more than enough. 
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Finally, the terms involving advection can be treated as: 


- (P • X,Pi-Ph) + (P ■ Xi,P ) - {P-XliPh) (±Xj) 

= ~{P ' (x~Xi),Pi - Ph ) - (P ■ XiiPi-Ph ) + {P ■ Xi,p) - (P-XhPh) 

= ~{P ■ (x~Xi),Pi - Ph ) + (P ■ Xi,P-Pi ) + ( Xi-XliPPh ) 

= ~{P ■ (x-Xi),Pi -Ph) + {P ■ Xi - P ■ Xi,P~Pi) + (Xi — Xii PPh - PPh) 
< Cph\\p! - Ph\\o(\\p ~ Pi\\o + lb-Philo + ft fe+ 1 |blk+i)- 


Inserting (99) (102) in (98) and using (54) and standard interpolation estimates we obtain (91) 


□ 


6 Numerical Experiments 


In this Section we will present some numerical experiments to validate the convergence results 
proven in the previous sections. We will test our method on the same problem and with the same 
meshes of |T3j , where we studied the Virtual Element Method for problem 0 in the primal form. 

Before presenting the numerical results we make a comment on the stabilization bilinear form 
in (44). For each element E £ Th we denote by xu f° r * = 1,2,.., Ne, the operator V^(E) —► R 
that to each Vh £ Vfi{E) associates the i-th local degree of freedom |I6|)-([I7|-(|T8|, ordered as 
follows: first the boundary d.o.f. (16), for i = 1,2, ...,iVj, and then the internal ones for 


i = Ng + 1,... ,Ne■ We assume that all the degrees of freedom are scaled in such a way that the 
associated dual basis l scales uniformly in the mesh size 


||0i|U°°(-E) — 1 Vi = 1,2,..., Ne- 


(103) 


With this notation, the most natural VEM stabilization S E (-,-) in (44) is given by (see nni) 


Ne 


s E (v -n° k v,w -n° k w) := \E\J2xi{v -n° k v) Xi(w- n°tu) 


(104) 


for all v,w £ V k (E). We now observe that, by definition of the L 2 projection operator 11°, and 
since both spaces Gk-i{E), Q k {E) appearing in (|T7|)-([T8|) are included in (P fe (E)) 2 , it is immediate 
to check that 

Xi (u-n°u)=0 Vu£V ft fc (E), i = N£ + l,...,N E . 


Therefore the contribution of the internal degrees of freedom in (104) vanishes, and we can equiv¬ 
alently use the shorter version 


-* V E 

S E (v-n° k v,w -n° k w) := |-E|^x*(v- n °v) Xi{w~n° k w). 

i—1 

In other words, the internal degrees of freedom do not need to be included in the stabilization 
procedure. 
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Figure 1: Lloyd-0 mesh 


Figure 2: Lloyd-100 mesh 


6.1 Exact Solution 

We will consider problem (JTJ) on the unit square with 

k(x, y) = (^_^y x 2 . b = {x,y), 7 = x 2 + y 3 , (105) 


and with right hand side and Dirichlet boundary conditions defined in such a way that the exact 
solution is 

p(x, y) = x 2 y + sin(27rx) sin(27ry) + 2. (106) 


The corresponding flux is given by 


u = —kS7p + bp. 


(107) 


We will show, in a loglog scale, the convergence curves of the error in L 2 between (p, u) and the 
solution (ph,Uh) given by the mixed Virtual Element Method (511. As the VEM flux Uh is not 
explicitly known inside the elements, we compare u with the L 2 —projection of onto (Pfc) 2 , that 
is, with 


6.2 Meshes 

For the convergence test we consider four sequences of meshes. 

The first sequence of meshes (labelled Lloyd-0) is a random Voronoi polygonal tessellation of 
the unit square in 25, 100, 400 and 1600 polygons. The second sequence (labelled Lloyd-100) is 
obtained starting from the previous one and performing 100 Lloyd iterations leading to a Centroidal 
Voronoi Tessellation (CVT) (see e.g. pH]). The 100-polygon mesh of each family is shown in Fig. [l] 
(Lloyd-0) and in Fig. [ 2 ] (Lloyd-100) respectively. 
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Figure 3: square mesh 


Figure 4: concave mesh 


The third sequence of meshes (labelled square) is simply a decomposition of the domain in 25, 
100, 400 and 1600 equal squares, while the fourth sequence (labelled concave) is obtained from 
the previous one by subdividing each small square into two non-convex (quite nasty) polygons. As 
before, the second meshes of the two sequences are shown in Fig. [3] and in Fig. [4] respectively. 

6.3 Convergence curves 

In Figs. [ 5 ] and [b] we report the relative error in L 2 for ph and Uh, respectively, for the four mesh 
sequences in the case k = 1. In Figs. [7] and [8] we show the same convergence results for k = 4. 

A closer inspection of the convergence curves for the L 2 error between p and ph shown in Figs. [ 5 ] 
and [7] reveals that the slope is slightly larger than expected for the coarsest meshes. This behavior 
can be explained in following way. The L 2 error \\p — ph\\o can be written as 

lb-Pfc|lo = l|P-P/Ho + lb/-Phllo (108) 

where we recall that on each element pj = II ® k p. As shown in Section [ 5 J there is a superconvergence 
of Ph to pp. 

||p/ “Philo < Ch k+2 . (109) 

Hence, as long as \\pj — Ph ||o is the dominant term in the error, we observe a slope of k + 2; when 
h becomes smaller, the term ||p — pi\\ 0 takes over and the slope becomes k + 1 as expected. This is 
clearly shown in Figs. [9| and |To| where p — pi and pi — p^ are plotted in the case of the lloyd-100 
meshes with k = 1 and k = 4, respectively. For the sake of clarity, on each curve we have reported 
its slope. 

We conclude that the Virtual Element Method behaves as expected and shows a remarkable 
stability with respect to the shape of the mesh polygons. 
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Figure 5: k = 1, relative L 2 error for ph 



Figure 6: k = 1, relative L 2 error for 


- * - Lloyd-0 



-O - Lloyd-100 



- □ - square 


* 

- V concave 




/ 

CD 


' //v 








vy 


*/& / 
<*'■ / 

5 

V 




mean diameter 



mean diameter 


Figure 7: k = 4, relative L 2 error for ph 


Figure 8: k = 4, relative L 2 error for Uh 
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Figure 9: k = 1, superconvergence 


Figure 10: k = 4, superconvergence 
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